### Table A 7 ###

rm(list = ls())
#install.packages("readstata13")
library("foreign")
library("readstata13")
library("dplyr")
library("sandwich")

data <- read.dta13("~/Dropbox/Replication_MVC/Datasets/datasets_analysis/panel_excombatientes.dta",
                   convert.factors = TRUE, generate.factors = FALSE,
                   encoding = "UTF-8", fromEncoding = NULL, convert.underscore = FALSE,
                   missing.type = FALSE, convert.dates = TRUE, replace.strl = TRUE,
                   add.rownames = FALSE, nonint.factors = TRUE, select.rows = NULL)

# Create member count
membercount <- data %>%
    arrange(year, wartimenetwork) %>%
    group_by(year, wartimenetwork) %>%
    summarize(members_total = n()) %>%
    ungroup()


# Merge the member count information back into the original dataset
data <- merge(data, membercount, by = c("year", "wartimenetwork"), all.x = TRUE)

vars <- c("codigoespejo", "origmun", "diffcaptures", "diffcaptures_f", "gold_shock_r", "year", "wartimenetwork", "oro_sin_titulo", "members_total", "Educationyears", "estratoeconomico", "region")
newdata <- data[vars]

# Separate by year
newdata2014 <- newdata[which(newdata$year=='2014'),]
newdata2015 <- newdata[which(newdata$year=='2015'),]
newdata2016 <- newdata[which(newdata$year=='2016'),]

# Dyad by year
dyad2014 <- merge(newdata2014, newdata2014, by = c("wartimenetwork"))
dyad2015 <- merge(newdata2015, newdata2015, by = c("wartimenetwork"))
dyad2016 <- merge(newdata2016, newdata2016, by = c("wartimenetwork"))

# Drop observation id.x=id.y
dyad2014n <- subset(dyad2014, codigoespejo.x!=codigoespejo.y)
dyad2015n <- subset(dyad2015, codigoespejo.x!=codigoespejo.y)
dyad2016n <- subset(dyad2016, codigoespejo.x!=codigoespejo.y)

#Production below the mean
dyad2014OST <- subset(dyad2014n, oro_sin_titulo.x<=39)
dyad2015OST <- subset(dyad2015n, oro_sin_titulo.x<=39)
dyad2016OST <- subset(dyad2016n, oro_sin_titulo.x<=39)

#Less than 100 in group
dyad2014NN <- subset(dyad2014n, members_total.x<=100)
dyad2015NN <- subset(dyad2015n, members_total.x<=100)
dyad2016NN <- subset(dyad2016n, members_total.x<=100)

# Dyad data
dyad <- rbind(dyad2014n, dyad2015n, dyad2016n)
dyadN <- subset(dyad, origmun.x!=origmun.y)

dyadOST <- rbind(dyad2014OST, dyad2015OST, dyad2016OST)
dyadOSTn <- subset(dyadOST, origmun.x!=origmun.y)

dyadNN <- rbind(dyad2014NN, dyad2015NN, dyad2016NN)
dyadNNn <- subset(dyadNN, origmun.x!=origmun.y)


# Analysis
# All
library("lfe")
Model1 <- felm(diffcaptures.x ~ gold_shock_r.y | year.x + Educationyears.x + estratoeconomico.x + region.x, data=dyadN)
Model2 <- felm(diffcaptures_f.x ~ gold_shock_r.y | year.x + Educationyears.x + estratoeconomico.x + region.x, data=dyadN)

Model3 <- felm(diffcaptures.x ~ gold_shock_r.y | year.x + Educationyears.x + estratoeconomico.x + region.x, data=dyadOSTn)
Model4 <- felm(diffcaptures_f.x ~ gold_shock_r.y | year.x + Educationyears.x + estratoeconomico.x + region.x, data=dyadOSTn)

Model5 <- felm(diffcaptures.x ~ gold_shock_r.y | year.x + Educationyears.x + estratoeconomico.x + region.x, data=dyadNN)
Model6 <- felm(diffcaptures_f.x ~ gold_shock_r.y | year.x + Educationyears.x + estratoeconomico.x + region.x, data=dyadNN)

library(stargazer)
stargazer(Model1, Model2, Model3, Model4, Model5, Model6,
          no.space = T,
          covariate.labels = c('Wartime Economic Shock'),
          keep.stat = c('n'))
